home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / TRED2.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  3KB  |  90 lines
  1. PROCEDURE tred2(VAR a: glnpnp; n: integer; VAR d,e: glnp);
  2. (* Programs using routine TRED2 must define the types
  3. TYPE
  4.    glnp = ARRAY [1..np] OF real;
  5.    glnpnp = ARRAY [1..np,1..np] OF real;
  6. where 'np by np' is the physical dimension of the matrix to be analyzed. *)
  7. VAR
  8.    l,k,j,i: integer;
  9.    scale,hh,h,g,f: real;
  10. FUNCTION sign(a,b: real): real;
  11.    BEGIN
  12.       IF (b < 0) THEN sign := -abs(a) ELSE sign := abs(a)
  13.    END;
  14. BEGIN
  15.    IF (n > 1) THEN BEGIN
  16.       FOR i := n DOWNTO 2 DO BEGIN
  17.          l := i-1;
  18.          h := 0.0;
  19.          scale := 0.0;
  20.          IF (l > 1) THEN BEGIN
  21.             FOR k := 1 TO l DO BEGIN
  22.                scale := scale+abs(a[i,k])
  23.             END;
  24.             IF (scale = 0.0) THEN BEGIN
  25.                e[i] := a[i,l]
  26.             END ELSE BEGIN
  27.                FOR k := 1 TO l DO BEGIN
  28.                   a[i,k] := a[i,k]/scale;
  29.                   h := h+sqr(a[i,k])
  30.                END;
  31.                f := a[i,l];
  32.                g := -sign(sqrt(h),f);
  33.                e[i] := scale*g;
  34.                h := h-f*g;
  35.                a[i,l] := f-g;
  36.                f := 0.0;
  37.                FOR j := 1 TO l DO BEGIN
  38.             (* Next statement can be omitted if eigenvectors not wanted *)
  39.                   a[j,i] := a[i,j]/h;
  40.                   g := 0.0;
  41.                   FOR k := 1 TO j DO BEGIN
  42.                      g := g+a[j,k]*a[i,k]
  43.                   END;
  44.                   IF (l > j) THEN FOR k := j+1 TO l DO g := g+a[k,j]*a[i,k];
  45.                   e[j] := g/h;
  46.                   f := f+e[j]*a[i,j]
  47.                END;
  48.                hh := f/(h+h);
  49.                FOR j := 1 TO l DO BEGIN
  50.                   f := a[i,j];
  51.                   g := e[j]-hh*f;
  52.                   e[j] := g;
  53.                   FOR k := 1 TO j DO a[j,k] := a[j,k]-f*e[k]-g*a[i,k]
  54.                END
  55.             END
  56.          END ELSE BEGIN
  57.             e[i] := a[i,l]
  58.          END;
  59.          d[i] := h
  60.       END
  61.    END;
  62.    (* Next statement can be omitted if eigenvectors not wanted *)
  63.    d[1] := 0.0;
  64.    e[1] := 0.0;
  65.    FOR i := 1 TO n DO BEGIN
  66.    (* Contents of this loop can be omitted if eigenvectors not wanted,
  67.       except for statement d[i] := a[i,i]; *)
  68.       l := i-1;
  69.       IF (d[i] <> 0.0) THEN BEGIN
  70.          FOR j := 1 TO l DO BEGIN
  71.             g := 0.0;
  72.             FOR k := 1 TO l DO BEGIN
  73.                g := g+a[i,k]*a[k,j]
  74.             END;
  75.             FOR k := 1 TO l DO BEGIN
  76.                a[k,j] := a[k,j]-g*a[k,i]
  77.             END
  78.          END
  79.       END;
  80.       d[i] := a[i,i];
  81.       a[i,i] := 1.0;
  82.       IF (l >= 1) THEN BEGIN
  83.          FOR j := 1 TO l DO BEGIN
  84.             a[i,j] := 0.0;
  85.             a[j,i] := 0.0
  86.          END
  87.       END
  88.    END
  89. END;
  90.